Anomalous force diffusion in nearly— ordered packings of frictionless discs 



cn 
O 
O 

>>: 



(N 



D.A. Headf 

^Department of Physics and Astronomy, Vrije Universiteit 1081 HV, Amsterdam, The Netherlands 

We derive analytic expressions for force propagation in packings of frictionless discs with a 
narrow distribution of disc sizes, by expanding to first order about the known ordered solution. The 
distribution of contact forces P(f) is found to be narrow at the upper surface, and broaden at a rate 
that varies with depth, being superdiffusive near the surface until crossing over to a subdiffusive 
regime near the fixed base. Furthermore, the response to an isolated load propagates along the 
edge of a 'cone,' as in the ordered case, but fluctuates under ensemble averaging by an amount that 
depends purely on height, not on the lateral position. Finally, we comment on ways in which the 
analytical framework presented here could be extended to a wider range of granular packings. 
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I. INTRODUCTION 



The mechanical properties of static bead packings re- 
mains a controversial and largely unanswered problem of 
granular media [1-4]. Experiments and numerical sim- 
ulations have thus far been unable to provide a defini- 
tive description of the propagation of stresses through 
rigid bead packings, in part because of the limited range 
of bead properties and system sizes typically considered. 
Theory is not yet at a sufficiently advanced stage to re- 
solve these deficiencies. First-principle approaches such 
as those initiated in [5-8] may eventually generate some 
form of 'complete' solution, but their significant complex- 
ity has thus far ruled out easily verifiable predictions. 
Analytically tractable theories have been devised, such 
as the g-model [9-11] and the force ray splitting 'dou- 
ble Y' model [3,12,13], but only after making ad hoc if 
intuitive assumptions regarding the local propagation of 
forces. 

An alternative approach to treating granular piles 
within an analytical framework is to restrict attention 
to a specific class of packings, but otherwise treat it as 
exactly as possible. Frictionless spheres are an ideal test 
case: the smooth surfaces cannot support tangential con- 
tact forces, removing indeterminacy problems that arise 
with frictional surfaces, and furthermore the spherical na- 
ture of the beads means that the remaining normal com- 
ponent must pass through the bead centre, so that the 
torque applied at each contact must vanish. Nonetheless 
disordered frictionless sphere packings are far from triv- 
ial: apart from mean field-type analysis [14], results ob- 
tained thus far have typically been numerical, either for 
small systems with arbitrary disorder [15,16], or much 
larger systems within the 'approximation of small dis- 
placements' of Roux et al. [17,18] (in this context, it is 
perhaps also relevant to mention the large-scale simula- 
tions of Breton et al. on geometrically ordered packings 
with disorder in the contact friction [19]). 

For a completely ordered lattice of frictionless discs, 
all of the forces propagate along the directions of col- 
inear bonds, and it is straightforward to write down all 



of the contact forces in response to a specified loading. 
In this paper, we expand about this ordered solution to 
first order in the deviations of bead radii from their mean 
value, and derive analytical expressions for the propaga- 
tion and distribution of forces for piles of arbitrary size. 
This we are able to eliminate statistical noise. Our find- 
ings are fully detailed in the subsequent sections, but 
in brief: the distribution of contact forces P(f) is nar- 
row near the surface, and broadens as one moves further 
down the pile. This broadening is initially superdiffusive 
(i.e. faster than the square root of distance from the 
surface), but crosses over to a subdiffusive regime lower 
down the pile and ceases to broaden precisely at the base. 
Turning to consider the Green's function problem of the 
response in contact forces to an isolated point load, we 
find that mean of this response is concentrated along a 
'cone' extending downwards from the point source along 
lines of contacts. Although the response inside the cone 
averages to zero, the variance in the fluctuations under 
ensemble averaging can be calculated, and is found to 
depend purely on the height, not (as might be expected) 
the distance from edge of the cone. 

This paper is arranged as follows. In Sec. II we spec- 
ify in detail the system to be considered. The statistical 
properties of the geometry of the packing under the speci- 
fied construction procedure and boundary conditions are 
derived in Sec. III. The force propagation from a sin- 
gle loaded bead within a particular packing is detailed 
in Sec. IV, which is then averaged under two types of 
loading to give the distribution of contact forces P(f) in 
Sec. V. In Sec. VI the response to an isolated load is 
found, including its statistical properties under ensem- 
ble averaging. Finally, in Sec. VII we summarise and 
suggest possible ways in which the analytical framework 
presented here might be extended. We remark that our 
approach is similar in spirit to Roux et al.'s 'assump- 
tion of small displacements,' although here were con- 
sider low-density packings for which the disorder can- 
not make or break contacts. In fact, it is essentially the 
small— polydispersity of the adaptive network algorithm 
of Tkachenko and Witten [16], treated analytically rather 
than numerically. 
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II. PRELIMINARIES 

Here we specify the class of granular packings under 
consideration in this paper, as well as the reasons for 
their analytically tractability. Firstly, we consider pack- 
ings that are marginally rigid or isostatic [14,20], i.e. the 
number of (scalar) constraints imposed by the conditions 
of force and torque balance on each bead equals the num- 
ber of degrees of freedom in the contact forces. Then the 
force network can be determined purely from the require- 
ments of local force and torque balance, without reference 
to the deformation properties of the beads. Also, since 
both the force and torque balance equations are linear, 
the full contact network can be found by simply summing 
the forces propagating from each loaded bead. 

Secondly, our packings are constructed from friction- 
less beads, so that only the normal components of the 
contact forces are non-zero. This means that the direc- 
tion of the contact forces can be uniquely determined 
from the geometry. Furthermore we use spherical beads, 
so that these contact forces point through the bead cen- 
tres and all torques trivially vanish. It is straightfor- 
ward to calculate the isostatic limit for such packings. 
If the mean number of contacts per bead is Z, this rep- 
resent Z/2 degrees of freedom in the (normal) contact 
forces, and d constraints imposed by force balance in a 
d-dimensional system. Thus the isostatic limit for fric- 
tionless spheres has Z = 2d. 

Thirdly, we suppose that the beads are almost 
monodisperse in size, i.e. the bead radii are narrowly 
distributed about some mean value. Then the contact 
topology will be the same as for a completely ordered sys- 
tem. Since Z = 2d everywhere, the load on any one bead 
can be uniquely decomposed into the d contact forces 
of the two lower beads (assuming there are no horizon- 
tal contacts). In this manner the forces due to external 
loads will propagate downwards through the system until 
reaching the fixed base [16,21]. 

Finally, the packing is formed by the sequential depo- 
sition of beads to the upper surface, and ensuring me- 
chanical stability between each addition. This exploits 
a happy coincidence of numbers: for Z = 2d packing, 
each incoming bead will rest on d others, which is also 
the number of contacts required to determine the posi- 
tion of the added bead. Thus the geometry of the pile 
can be found without needing to know the incoming bead 
trajectories. 



III. GEOMETRY OF THE PACKING 

In this section, the position vectors of the beads for a 
particular set of bead radii is derived for the d = 2 case 
of frictionless discs. An assumption underlying the anal- 
ysis throughout this paper is that all of the contact forces 
are compressive. This is crucial, since granular particles 



are typically regarded as non-cohesive and therefore can- 
not support tensile contacts. It is clear that such a state 
must exist if the disorder is sufficiently small, since all 
contact forces are compressive in the corresponding or- 
dered limit [22]. The range of validity of this assumption 
will be discussed in Sec. V after the force distribution has 
been found. 

Before considering the near-crystalline packing, it is 
useful to describe the completely ordered case, both as 
reference and to fix the notation. See Fig. 1. Without 
disorder, each bead has radius r and is located at a lat- 
tice site (i,j), where i and j refer to the horizontal and 
vertically-upwards lattice directions, respectively. The 
centre of the bead at is denoted by the position 

vector Xjj , where 

x ij = r U + + r(j - i)n~ (1) 

Here, are the unit base lattice vectors in the upwards- 
right (n+) and upwards-left (n~) directions. It is 
straightforward to show that n ± = (±n x ,n y ) with 
n x = s/Ar, where s is the horizontal distance between 
bead centres. Each bead makes 2d = 4 contacts with 
beads aligned along the lattice diagonals, i.e. bead (i, j) 
touches beads (i ± l,j ± 1). Thus only sites with i + j 
even are occupied. The requirement that beads cannot 
overlap in either the horizontal or vertical directions fixes 
s to be in the range (2r,2\/3r). 

Disorder is introduced into the packing by assigning 
each bead a radius r + Srij prior to deposition onto lat- 
tice site (i, j), where the Sr^ are random variables drawn 
from a distribution with zero mean and variance a\ r . 
When a bead is added to the growing surface, its cen- 
tre Xy + (5xy must be chosen to ensure the distance to 
both supporting beads is equal to the sum of the corre- 
sponding radii; to first order in the Sr^ , this gives 

1/0 %L \ 

+ (Sr i:j + (5ri_ij_i) b + + (Sr^ + 5r l+1 b (2) 

where b* are the reciprocal lattice vectors b* = 
(±l/2n x , l/2n y ), i.e. vectors b^ 1 such that n 1 * 1 • b ± = 1 
and n ± • b^ = 0. 

We now fix the lower boundary condition to be that the 
centres of all beads in the base j = have the same height 
and are uniformly spread horizontally: x^o = (2irn x ,0), 
just as in the ordered case. Then (2) can be extrapolated 
to j = via an induction process to give the perturbed 
bead positions purely in terms of the Sr^ , 

5xij = < 5n-jfi + Snj + 2 ^2 Sr i-U-m),m \ b + 

I m=l J 

+ |^+j,o + Snj +2^2 <fr*+0'-™)."» j b_ ( 3 ) 
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Thus the position of the bead at depends on 5rij, 

5ri±i t j—i, 5r>i±2.j-2 and so on, i.e. only on the Srij ly- 
ing along the two diagonals stretching from to the 
base, hereafter referred to as the 'cone.' Note that (3) ig- 
nores any horizontal boundaries, and is therefore valid for 
infinitely wide 'slabs' and systems with periodic bound- 
aries at i = ±W with W > M, where j = M denotes the 
maximum height of the packing. 

We can now describe the manner in which the geo- 
metric disorder propagates upwards from the base. For 
j > 1, the expression (3) is dominated by the two sums. 
Thus, to good approximation, the coefficients of both 
the b + and b~ components of 5xy are the sum of j 
independent random variables 2Sr k[ . Each such term 
has a variance 4cr| r , and so the combined variance is 
~ 4:j<Jg r . Therefore the x and y-components of <5xy will 
both broaden like O^^crgr), i- e - the geometric disorder 
diffuses upwards in a normal (Brownian) manner. This 
contrasts with the non-Brownian diffusion of forces de- 
scribed in Sec. V. 

For frictionless discs, the contact force between any 
two contacting beads is parallel to the line connecting 
the bead centres, which can be derived using (3). For 
instance, if the direction from the centre of bead to 
(i + l,j + 1) is denoted n + + 5nt, then 

2rSh+ = (6nj + Sr i+lj+1 )(b + + b n+) 

+ b j(5ri + j + 2,o - Sr l+j . 



+ 2 y^[^ r i+(j-m)+2 ; n J u 
rn—1 



(4) 



This obeys <5n^ • n+ = 0, ensuring that n + + 5n^ is a 
unit vector to first order in Sn^. A similar expression 
applies for contacts from to (i — l,j + 1), which 

have direction n~ + Shf- with 

2rSh~ j = (Snj + 5rj_ij+i)(b + + b n~) 

+ b+|(5r i _ i _ 2 ,o - Sn-jfi 

+ 2 X! [^ r i-(i-m)-2,m - ^i-(j-m),m] j 



m—1 



IV. 



PROPAGATION OF FORCES FROM A 
SINGLE LOADED BEAD 



Consider the force propagating from a vertical load 
F ioad = ( Q) _yioad) applied t0 a single bead ;) If the 

packing were ordered, i.e. Sr^ = everywhere, F load 



would first split equally into both of the bead's support- 
ing contacts, and then propagate unaltered along colinear 
bonds to the fixed base. With disorder Sr^ ^ 0, there 
are two differences: not only is F load shared unequally 
by the two contacts (which now generally have different 
directions), but also the propagation to the base is no 
longer linear but multiply branched. A qualitative de- 
scription of the following analysis can be found in [23] 
for close-packed 3D crystals. 

Let denote the magnitude of the contact force from 
(i,j) to (i ± 1, j + 1), using the convention that com- 
pressive contacts have positive f^. The force balance 
equation for the loaded bead (k, I) is then 



riload 



+ (n + +Sn+_ U -i)fLil-i 
+ (»~ +tnf+ii-i)fk+u-i 



(6) 

By writing F load as -(/ load /2n a )(n+ + n"), the support- 
ing forces can be evaluated to first order in the Sn^, 



fload 



n • Shf ,, , - n+ • Sh 



f+ _ /""" J i , ' ""A- :/■ I 

Jk-ll-1 — ~ ^ 1 + 



k+ll-1 



2n, 



1 - , 



,- _ / load f <M+ ■ 5hf +ll _ 1 - n ■ <Sn+_n-i 1 
Jk+U-i- 2n v Y + 1-02 j 



(7) 



using (f> = n+ • n~. The range of allowed s given in 
Sec. HI means that — \ < (j) < ^, so that 1 — (j) 2 > | and 
the denominators in (7) never diverge. 

Once the two contact forces f k ± Tll _ 1 are known, it is 
necessary to calculate how they propagate to the fixed 
base. This is done by the same force balance equation as 
(6) but with F load replaced by the corresponding incom- 
ing contact force. For instance, the force /+ coming into 
a bead (r, s) along the — (n + + 5n+ ) direction is balanced 
by the two supporting contacts f^ fl which to order 
0(5^) are 

^f i = - T 1 ^n--(Sh+_ ls - 1 -Sh+) (8) 

The expressions for the contributions from — /^(n - + 
Shf s ) follow from symmetry. Thus the contribution of 
each contact force of the loaded bead /^ pl l _ 1 to any 
(5) contact in the system, if any, can be found by sum- 
ming over all downward-propagating paths from (k, I) to 
(i,j), and applying (8) for each step. This process is 
simplified by the observation that the "branching" force 
fr+is-i m (8) is 0(Sh ± ) to leading order, and thus only 
paths with or 1 branch point need to be considered 
within this first order calculation. 

Without loss of generality we now consider the prop- 
agation of the fk-ii-i forces; the equivalent f k ~ +li _ 1 re- 
sults follow from symmetry. Suppose a contact ff- lower 
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down the packing obeys k — i = I — j , so that it is con- 
nected to along a straight path with no branches. 
Then (8) can be applied iteratively to give 



f + - 



fk-ll-l 



= 1 + 



1 



r n-.(Jn+-«yn+_ 1J _ 1 ) (9) 



Only the <5n + at each end of the path contribute: the 
intermediate terms cancel, to first order. 

Contacts f~j with the opposite orientation may still re- 
ceive a contribution from f^-u-i ^ a P a *h with precisely 
one branch point connects the two bonds. This happens 
if n = i [(/ — j) — (k — i)] obeys 1 < n < I — j, in which 
case the branch point is at bead (i — n,j + n). As the 
branch is 0(5^) to leading order, only the zero'th order 
contributions from the remaining propagation terms are 
required, giving 



1 



fk-ll-l 



(10) 



All other /± receive zero contribution. Thus the con- 
tribution to any contact force f^j due to the load F load 
applied to (k, I) can be found by combining the expres- 
sions (7), (9) and (10). 



V. DISTRIBUTION OF FORCES P(F) 



When a specified load is applied to the pile, the total 
contribution to any contact force in the packing can 
be found by applying the results of Sec. IV to each loaded 
bead, and summing. The distribution of contact forces 
P(f) can then be found by either averaging over all con- 
tacts with the same height, or by ensemble averaging over 
different packings {Srij}; for the loads considered here, 
which do not vary with horizontal distance, the final re- 
sult is the same. Two alternative loadings are analysed: 
a surface load F load = (0, - f load ) applied to all beads 
at j '■ — M with all other beads weightless; and a bulk 
load in which each bead has a weight (0, —vug). Only the 
first case will be treated in detail here; the results of the 
qualitatively-similar bulk loading will be simply stated 
at the end. 

Each contact ff- will have contributions from the sur- 
face bead (i + M — j, M) according to the unidirectional 
propagator (9), plus another M—j—1 contributions from 
beads in the range (i-M+j + 2,M) to (i + M-j-2,M) 
which branch once according to (10). Including the initial 
splitting of F load at the surface (7), the final expression 
for f± is 



F 
j i 



1 + ' ^ ~ ^ ' 5A i+M-j+l,M-l) 

M-j-1 

n—1 

In principle, it is possible to write this expression in terms 
of the Sr^; however, this soon becomes messy due to 
the proliferation of terms. Instead we now assume that 
M — j ;» 1, so that (i, j) is 'sufficiently' below the sur- 
face. In this limit, the second term on the right hand 
side of (11) is much smaller than the sum and can there- 
fore be dropped (recall that Siifj ~ C(j 1//2 )). Further- 
more, the sum itself will be dominated by terms with 
n = 0(M — j), so we need only consider 5n~ s with 
s = 0(M — j). This means that the sum terms in the ex- 
pressions for Sn^ (4,5) will dominate and the remaining 
terms can be dropped. Combining these approximations, 
f£ can be written purely in terms of combinations of Sr rs 
like 2<5rv 



Sr. 



r-2,. 



8r 



r+2,s 



A r . 



f + - 



^ M-j j+n 

f™ d /2n y r(l - <}> 2 ) S S A i-U-m),m 

■> i a n=l m=l 

M M-j 

m=ln=max(l,m-i) 



r{l-4> 2 ) 



M 

E 



M — j : m < j 
M — m : m > j 



A, 



i — (j~m),m 



(12) 



Since the mean of each A rs is zero, we can immedi- 
ately see that the mean force (/) = f load /2n y for all 
heights. Also, since the A i _(j_ m \ m are independent for 
each value of m, then /rt is the sum of a number of inde- 
pendent random variables, each of which has a weighting 
of no more than 0(1/ M) of the total. Hence the central 
limit theorem applies, and, assuming that the variance of 
bead radii aj r is finite, the distribution of contact forces 
must be Gaussian. The variance of forces C/(j) depends 
on the height j and can be found by standard techniques, 
given the weightings of the A i _^_ m - ) m shown in (??). 
Noting that the variance of A rs is 6(T^ r , we finally find 
that 



(f) 2 



2M 3 



1 



h2\2 



(l-z) 2 (2z + l)(^) 2 (13) 



/ load /2n y 



where z — j/M. A graphical representation of this so- 
lution is given in Fig. 2. This should be compared with 
the Brownian expectation <r| ~ 1 — z, as seen in e.g. the 
g-model [10,11]. Not only is (13) superdiffusive near the 
surface, cr 2 ~ (1 — z) 2 , the prefactor depends on M, i.e. 
the depth of the pile is 'felt' at the surface, no matter how 
large M is. This is not unexpected, since the degree of 
geometric disorder at the surface always depends on M in 
these near-crystalline packings; for more disordered sys- 
tems, it would most likely saturate at some finite height 
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and the M-dependence near the surface would vanish. 
Also note that the expression (13) becomes subdiffusive 
nearer the base, and ceases to broaden precisely at j = 0. 
This is because the geometric disorder vanishes at the 
base, so that the forces no longer diffuse. To check the 
assumptions leading to (13), we compare it to data ob- 
tained from the adaptive network algorithm [16] in Fig. 3. 
The agreement is good, and appears to improve as M is 
increased. 

For a bulk loaded system, the above calculations can 
be repeated in a similar manner. The main difference 
is that the mean force now varies with depth, (f) z = 
Mmg(l — z)/2n y . Nonetheless, the variance normalised 
by the mean force for each depth takes a similar form, 



3M 3 



^(l-z) 2 (4z + l)(^) 2 (14) 



</>2 ~ 10(1-02) 

The right hand side of this expression reaches a maximum 
at finite height j* = M/6, and then narrows slightly to- 
wards to the base. This is because the diffusion of forces 
decreases near the base, as mentioned above, but the 
mean force (to which the distribution is normalised) in- 
creases linearly with 1 — z, which acts to narrow the dis- 
tribution. If the distribution was not normalised to (f) z , 
then it would broaden monotonically and again cease to 
broaden just at the base, for the same reasons as the 
surface loading case. 



A. Breakdown of the compressive bond assumption 

We can now discuss the range of validity of the assump- 
tions made earlier. The magnitude of the displaced bead 
centres 5xij typically vary as 0{M 1 / 2 crs r ), which will be 
similar to the bead radii r, and hence violate the as- 
sumed contact topology, for piles of height 0[(asr l r )~ 2 ] 
and greater (this assumes that the gaps between beads 
in the ordered packing are also 0(r)). However, the as- 
sumption of non-tensile contacts will be violated much 
sooner. For both surface and bulk loading, 07 averaged 
over the system has the same general form 



0~6r 



(/> r 



(15) 



Thus a finite fraction of forces will become negative for 
packings of height M* = 0[(as r /r)~ 2 / 3 ], and the above 
analysis only holds for M <C M*. Note that this is signif- 
icantly smaller than if the forces diffusion in a Brownian 
manner ~ M 1 ' 2 . Strictly speaking, we should really re- 
quire that a vanishing number of contact forces are nega- 
tive rather than a vanishing fraction, which would give a 
much lower M* . However, this is most likely too restric- 
tive, unless it could be shown that e.g. a single negative 
force initiates a cascade of rearrangements that alters the 
forces throughout a finite fraction of the system. 



VI. RESPONSE GREEN'S FUNCTION 

The various phenomenological theories for granular 
stress propagation differ most markedly in the response 
Green's function to an isolated load [1-3,24]. This can 
be calculated for our system as follows: apply a load 
(0, — J load ) to bead (k, I), and measure the vertical com- 
ponent p of the induced force at a contact / 4 • . This 
can be found for a given packing {Srij } using the results 
given above. Now ensemble average over different pack- 
ing realisations, i.e. uncorrelated sets of {Sr^} with the 
same variance a 2 r . Since this is a first-order calculation, 
the mean response will be as the ordered case, so that 
(p) ensemble will be zero except along the cone propagating 
downwards from (k,l). However, although the response 
outside this cone is strictly zero, the response inside the 
cone only vanishes after averaging: it will be positive or 
negative for each particular packing. The variance of the 
distribution of p between different realisations can be cal- 
culated using a similar procedure to that given in Sec. V, 



ensemble (P) 
(/ 



load ^2 



2(1-02): 



0~,h 



(16) 



A schematic representation of this solution is given in 
Fig. 4. Interestingly, this does not decay away from 
the edges of the cone, as predicted by phenomenologi- 
cal approaches, and in fact is completely independent of 
the horizontal position, depending only on the combined 
heights I and j. 

The counterpart to ensemble averaging is coarse grain- 
ing, i.e. applying a uniform pressure to a group of N 
beads, and jointly taking N and j to infinity in a suit- 
able manner to give a continuum result. However, for our 
packings this will just give a zero response everywhere ex- 
cept along the edge of the cone. This is because we have 
only expanded around the ordered solution to first order, 
which will always give the zero-th order, ordered result 
after any form of averaging. Extending these calcula- 
tions to second order will introduce quadratic terms that 
do not vanish under averaging, and thus would allow a 
direct comparison between different types of averaging. 



VII. SUMMARY AND DISCUSSION 

In summary, we have obtained analytical expressions 
for contact forces in isostatic packings of frictionless discs 
in the nearly ordered limit. These explicitly demon- 
strate the anomalous broadening of the distribution of 
forces P(f) with increasing depth, and the unusual re- 
sponse Green's function, whose magnitude of fluctua- 
tions under ensemble averaging depends only on the ver- 
tical coordinate. We see no reason in principle why 
these findings could not be tested experimentally, as long 
as a sufficiently monodisperse sample of smooth beads 
could be found. Indeed, this may be the main prob- 
lem: glass bead experiments with a low-polydispersity 
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of sizes placed in close-packed crystalline configurations 
have shown that P(f) already has a broad exponential 
P(f) at the base [23,25], as opposed to the Gaussian form 
predicted here. This suggests that even these carefully 
controlled experiments are already too disordered to see 
the first-order results predicted here. 

It is hoped that the calculations presented in this 
paper may lay the foundations for a more complete, 
analytically-tractable theory of bead packings that 
would complement current approaches. Thus it is worth 
discussing in what ways these calculations may be ex- 
tended to more realistic situations. There should be lit- 
tle problem in going to 3 dimensions, or for considering 
loads at an angle to the vertical, which would allow the 
piles to be sheared. Perhaps more interesting would be 
to expand about the ordered solution to second order 
in the Srij, rather than just to first order as considered 
here. This would allow the mean force to deviate from 
the crystalline solution for the first time, and hence al- 
low the nature of force propagation in granular packings 
to be addressed from within an analytical framework; 
for instance, it would be possible to see if proposed re- 
lationships between the components of the stress tensor 
are valid [1,14]. It would also show if reaching the con- 
tinuum limit by ensemble averaging differs from coarse 
graining, as already discussed in Sec. VI. 

However, extending this analysis to frictional beads is 
likely to be more troublesome. The marginal rigidity 
state for beads with infinite friction has Z = 3 in 2 di- 
mensions, so that sequentially deposited beads will make 
either 1 or 2 contacts on the surface - in the case of 1 con- 
tact, it is not possible to determine the rest position of 
the incoming bead. Thus groups of beads must be added 
simultaneously, so that the final geometry will depend on 
incoming bead trajectories and their material properties. 
Also, with Z — 3 it is not possible to uniquely decompose 
the load applied to any given bead along paths connect- 
ing it to the base, for either the normal or tangential con- 
tact forces. Instead loops, and also paths leading to the 
surface, will inevitably arise. This non-unique decompo- 
sition will also be a problem for frictionless non-spherical 
beads, for which Z = 6 in 2 dimensions. It is not clear 
how these problems may be surmounted without resort- 
ing to ad hoc assumptions or simplifications. 
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FIG. 1. An ordered packing of discs of radius r and hori- 
zontal separation s (thick black circles). The lattice indices 
and the lattice base vectors are also shown. A par- 
ticular example of a disordered packing, in which the beads 
have radii r + Srij , is superimposed in grey. 
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(a) (b) 
FIG. 2. Graphical representation of (a) the broadening of 
the geometric disorder with height j, and (b) the broadening 
of the force distribution with distance from the surface, as 
given by the surface-loaded solution (13). 
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FIG. 3. Comparison of the theoretical prediction for force 
diffusion for a surface-loaded system (13) with corresponding 
simulation data. The normalised standard deviation of forces 
Of /(f) has been scaled by M 3//2 so that, according to (13), 
the curve will be M-independent when plotted against the 
relative height z = j/M. The two pile heights used in the 
simulations are given in the key, and appear to confirm this 
prediction. The parameters for both theory and simulations 
are (a Sr /r) 2 = lCT 6 / 12 , r = 1 and s = 3. 




(a) (b) 
FIG. 4. Schematic representation of the Green's response 
solution (16) to a vertical load F load . (a) The mean re- 
sponse is identical to the corresponding ordered system, in 
which forces propagate along lines of colinear bonds forming 
a 'cone.' (b) The magnitude of the fluctuations under ensem- 
ble-averaging, with darker grey corresponding to a greater 
variance. The variance is zero outside the cone, depends only 
on height inside the cone, and is twice as large near the apex as 
near the base (fluctuations on the cone itself are not shown) . 
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